#include "RayCasting.h"
#include "Ray.h"
#include "png.h"
#include <math.h>

int normalizeVectoriR(float * e)
{
        float magnitude = sqrtf(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);

        e[0] /= magnitude;
        e[1] /= magnitude;
        e[2] /= magnitude;

        return 0;
}



float InterpolacionTri(volume_t * volume, float * p)
{
#if 0
        int x1 = (int) p[0];
        int y1 = (int) p[1];
        int z1 = (int) p[2]; 
        int x2 = x1+1;
        int y2 = y1+1;
        int z2 = z1+1;

        return  volume->data[x1][y1][z1] * (1.0-p[0]) * (1.0-p[1]) * (1.0-p[2]) +       \
                volume->data[x2][y1][z1] * p[0]       * (1.0-p[1]) * (1.0-p[2]) +       \
                volume->data[x1][y2][z1] * (1.0-p[0]) * p[1]       * (1.0-p[2]) +       \
                volume->data[x1][y1][z2] * (1.0-p[0]) * (1.0-p[1]) * p[2]       +       \
                volume->data[x2][y1][z2] * p[0]       * (1.0-p[1]) * p[2]       +       \
                volume->data[x1][y2][z2] * (1.0-p[0]) * p[1]       * p[2]       +       \
                volume->data[x2][y2][z1] * p[0]       * p[1]       * (1.0-p[2]) +       \
                volume->data[x2][y2][z2] * p[0]       * p[1]       * p[2];
#endif
#if 1
        int x0 = floorf(p[0]);
        int y0 = floorf(p[1]);
        int z0 = floorf(p[2]);
	if (x0 < 0.0) x0 = 0.0; if (y0 < 0.0) y0=0.0; if(z0 <= 0.0)z0=0.0;
	if (x0 >= volume->dim[0]) x0 = volume->dim[0]-2; if (y0 >= volume->dim[1]) y0=volume->dim[1]-2; if(z0 >= volume->dim[2])z0=volume->dim[2]-2;
        int x1 = x0 + 1 < volume->dim[0] ? x0 + 1 : x0;
        int y1 = y0 + 1 < volume->dim[1] ? y0 + 1 : y0;
        int z1 = z0 + 1 < volume->dim[2] ? z0 + 1 : z0;

//      fprintf(stdout, "(%f,%f,%f) (%d,%d,%d)  (%d,%d,%d)\n",p[0],p[1],p[2],x0,y0,z0,x1,y1,z1);
/*
        float p000 = volume->data[x0][y0][z0];
        float p001 = volume->data[x0][y0][z1];
        float p010 = volume->data[x0][y1][z0];
        float p011 = volume->data[x0][y1][z1];
        float p100 = volume->data[x1][y0][z0];
        float p101 = volume->data[x1][y0][z1];
        float p110 = volume->data[x1][y1][z0];
        float p111 = volume->data[x1][y1][z1];
*/
        float p000 = volume->data[x0][y0][z1];
        float p001 = volume->data[x0][y1][z1];
        float p010 = volume->data[x0][y0][z0];
        float p011 = volume->data[x0][y1][z0];
        float p100 = volume->data[x1][y0][z1];
        float p101 = volume->data[x1][y1][z1];
        float p110 = volume->data[x1][y0][z0];
        float p111 = volume->data[x1][y1][z0];

	float pi[3];
	pi[0] = p[0]-x0;
	pi[1] = p[1]-y0;
	pi[2] = p[2]-z0;

/*	
fprintf(stdout,"(%f,%f,%f) --> (%d,%d,%d)=%f (%d,%d,%d)=%f (%d,%d,%d)=%f (%d,%d,%d)=%f (%d,%d,%d)=%f (%d,%d,%d)=%f (%d,%d,%d)=%f (%d,%d,%d)=%f\n",		p[0],p[1],p[2],				\
																			x0,y0,z0,volume->data[x0][y0][z0],	\
																			x0,y0,z1,volume->data[x0][y0][z1],	\
																			x0,y1,z0,volume->data[x0][y1][z0],	\
																			x0,y1,z1,volume->data[x0][y1][z1],	\
																			x1,y0,z0,volume->data[x1][y0][z0],	\
																			x1,y0,z1,volume->data[x1][y0][z1],	\
																			x1,y1,z0,volume->data[x1][y1][z0],	\
																			x1,y1,z1,volume->data[x1][y1][z1]);

fprintf(stdout,		"V000(1-x)(1-y)(1-z)\t\t=%f*%f*%f*%f\t\t=%f\n"	
			"V001(1-x)(1-y)z\t\t=%f*%f*%f*%f\t\t=%f\n"	
			"V010(1-x)y(1-z)\t\t=%f*%f*%f*%f\t\t=%f\n"	
			"V011(1-x)yz\t\t=%f*%f*%f*%f\t\t=%f\n"	
			"V100x(1-y)(1-z)\t\t=%f*%f*%f*%f\t\t=%f\n"	
			"V101x(1-y)z\t\t=%f*%f*%f*%f\t\t=%f\n"	
			"V110xy(1-z)\t\t=%f*%f*%f*%f\t\t=%f\n"	
			"V111xyz\t\t=%f*%f*%f*%f\t\t=%f\n",
			p000,(1.0-p[0]),(1.0-p[1]),(1.0-p[2]),p000 * (1.0-p[0]) * (1.0-p[1]) * (1.0-p[2]),\
			p001,(1.0-p[0]),(1.0-p[1]),(p[2]),p001 * (1.0-p[0]) * (1.0-p[1]) * p[2],\
			p010,(1.0-p[0]),(p[1]),(1.0-p[2]),p010 * (1.0-p[0]) * p[1]       * (1.0-p[2]),\
			p011,(1.0-p[0]),(p[1]),(p[2]),p011 * (1.0-p[0]) * p[1]       * p[2],\
			p100,(p[0]),(1.0-p[1]),(1.0-p[2]),p100 * p[0]       * (1.0-p[1]) * (1.0-p[2]),\
			p101,(p[0]),(1.0-p[1]),(p[2]),p101 * p[0]       * (1.0-p[1]) * p[2],\
			p110,(p[0]),(p[1]),(1.0-p[2]),p110 * p[0]       * p[1]       * (1.0-p[2]),\
			p111,(p[0]),(p[1]),(p[2]),p111 * p[0]       * p[1]       * p[2]);
*/

        return  p000 * (1.0-pi[0]) * (1.0-pi[1]) * (1.0-pi[2]) +       \
                p001 * (1.0-pi[0]) * (1.0-pi[1]) * pi[2]       +       \
                p010 * (1.0-pi[0]) * pi[1]       * (1.0-pi[2]) +       \
                p011 * (1.0-pi[0]) * pi[1]       * pi[2]       +       \
                p100 * pi[0]       * (1.0-pi[1]) * (1.0-pi[2]) +       \
                p101 * pi[0]       * (1.0-pi[1]) * pi[2]       +       \
                p110 * pi[0]       * pi[1]       * (1.0-pi[2]) +       \
                p111 * pi[0]       * pi[1]       * pi[2];
#endif


}


int Intersection(volume_t volume, float * org, float * dir, float * startP, float * endP)
{
	float dx 	= -org[0]/dir[0];
	float dy 	= -org[1]/dir[1];
	float dz 	= -org[2]/dir[2];
	float xdx	= (volume.dim[0] - org[0])/dir[0];
	float ydy	= (volume.dim[1] - org[1])/dir[1];
	float zdz	= (volume.dim[2] - org[2])/dir[2];

	float p1 = 0;
	float p2 = 0;
	float p3 = 0;
	float p4 = 0;
	float p5 = 0;
	float p6 = 0;

	if ((dir[0] == 0) && (dir[1] == 0) && (dir[2] == 0))
	{
		return 0;		
	}
	else if ((dir[0] == 0) && (dir[1] == 0) && (dir[2] != 0))
	{
		return 0;
	}
	else if ((dir[0] == 0) && (dir[1] != 0) && (dir[2] == 0))
	{
		return 0;
	}
	else if ((dir[0] == 0) && (dir[1] != 0) && (dir[2] != 0))
	{
		return 0;
	}
	else if ((dir[0] != 0) && (dir[1] == 0) && (dir[2] == 0))
	{
		return 0;
	}
	else if ((dir[0] != 0) && (dir[1] == 0) && (dir[2] != 0))
	{
		return 0;
	}
	else if ((dir[0] != 0) && (dir[1] != 0) && (dir[2] == 0))
	{
		return 0;
	}
	else
	{
		/* 
		 * Plano 
		 * (0,0,0)-(x,0,0)
		 * (x,0,0)-(x,0,z)
		 * (x,0,z)-(0,0,z)
		 * (0,0,z)-(0,0,0)
		 */
		if ((dy >= dx) && (dy <= xdx) && (dy >= dz) && (dy <= zdz))
			p1 = dy;

		/* 
		 * Plano 
		 * (0,0,0)-(x,0,0)
		 * (x,0,0)-(x,y,0)
		 * (x,y,0)-(0,y,0)
		 * (0,y,0)-(0,0,0)
		 */
		if ((dz >= dx) && (dz <= xdx) && (dz >= dy) && (dz <= ydy))
			p2 = dz;

		/* 
		 * Plano 
		 * (0,0,0)-(0,y,0)
		 * (0,y,0)-(0,y,z)
		 * (0,y,z)-(0,0,z)
		 * (0,0,z)-(0,0,0)
		 */
		if ((dx >= dy) && (dx <= ydy) && (dx >= dz) && (dx <= zdz))
			p3 = dx;

		/* 
		 * Plano 
		 * (x,0,0)-(x,y,0)
		 * (x,y,0)-(x,y,z)
		 * (x,y,z)-(x,0,z)
		 * (x,0,z)-(x,0,0)
		 */
		if ((xdx >= dy ) && (xdx <= ydy) && (xdx >= dz) && (xdx <= zdz))
			p4 = xdx;

		/* 
		 * Plano 
		 * (0,y,0)-(x,y,0)
		 * (x,y,0)-(x,y,z)
		 * (x,y,z)-(0,y,z)
		 * (0,y,z)-(0,0,z)
		 */
		if ((ydy >= dx) && (ydy <= xdx) && (ydy >= dz) && (ydy <= zdz))
			p5 = ydy;

		/* 
		 * Plano 
		 * (0,0,z)-(x,0,z)
		 * (x,0,z)-(x,y,z)
		 * (x,y,z)-(0,y,z)
		 * (0,y,z)-(0,0,z)
		 */
		if ((zdz >= dx) && (zdz <= xdx) && (zdz >= dy) && (zdz <= ydy))
			p6 = zdz;
	}

	fprintf(stdout,"%f %f %f %f %f %f",p1,p2,p3,p4,p5,p6);

	

	return 0;
}

int reducingVector(float * e)
{
        float magnitude = sqrtf(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);

        e[0] /= magnitude;
        e[1] /= magnitude;
        e[2] /= magnitude;

        return 0;
}

float Intersection2(volume_t volume, float * org, float * dir, float * ed, unsigned char * s)
{
	int i = 0;
	int p[3];
	float pi[3];
	float put = 0;
	int   m = 0;
	int primera = 1; 
	float ant, sig;
	float iso = 0.297761;

	for(i=0; i<1000000; i++)
	{
		
		p[0]=(int)org[0] + i*dir[0];
		p[1]=(int)org[1] + i*dir[1];
		p[2]=(int)org[2] + i*dir[2];
		
		pi[0]=org[0] + i*dir[0];
		pi[1]=org[1] + i*dir[1];
		pi[2]=org[2] + i*dir[2];
		
		if ((p[0]>=0)&&(p[0]<volume.dim[0])&&(p[1]>=0)&&(p[1]<volume.dim[1])&&(p[2]>=0)&&(p[2]<volume.dim[2])){
//	fprintf(stdout,"%f %f %f %d %d %d --->\n",pi[0],pi[1],pi[2],p[0],p[1],p[2]);
			put += p[1];//volume.data[p[0]][p[1]][p[2]];
			if (primera)
			{
				ant = InterpolacionTri(&volume,pi);//volume.data[p[0]][p[1]][p[2]];
				primera = 0;	
			}
			else
			{
				sig = InterpolacionTri(&volume,pi);//volume.data[p[0]][p[1]][p[2]];
				if (( ((iso-ant)<0.0) && ((iso-sig)<0.0)) || ( ((iso-ant)>0.0) && ((iso-sig)>0.0)))
				{
					ant = sig;
				}
				else
				{
//			fprintf(stdout,"tocado (%f,%f,%f)=%f \n",pi[0],pi[1],pi[2], InterpolacionTri(&volume,pi));
					*ed=sig;
//					normalizeVectoriR(pi);
					s[0]=(pi[0]/volume.dim[0])*255;
					s[1]=(pi[1]/volume.dim[1])*255;
					s[2]=(pi[2]/volume.dim[2])*255;
					return p[1];
				}
			}
			m++;
		}
	}

	if (!(i<1000000))
		return 0.0;
	
}

int IntersectionPlane(volume_t volume, float * org, float * dir, float * startP, float * endP)
{
	
	float asd;
	/* plano 1 */
return 0;
	

}

int rayCasting(volume_t volume, Camera_t camera, Screen_t screen)
{

	int i = 0;
	int j = 0;
	int k = 0;

	float ds;
	float de;

//	Ray_t ray;
	float ph[3];
	float pw[3];
	float sr[3];
	float dr[3];
	float ih  = screen.h/screen.H;
	float iw  = screen.w/screen.W;
	sr[0] = camera.position[0];
	sr[1] = camera.position[1];
	sr[2] = camera.position[2];

	ph[0] = screen.distance * (tan(screen.fov_H*(M_PI/180.0)) * camera.up[0]);
	ph[1] = screen.distance * (tan(screen.fov_H*(M_PI/180.0)) * camera.up[1]);
	ph[2] = screen.distance * (tan(screen.fov_H*(M_PI/180.0)) * camera.up[2]);
	pw[0] = screen.distance * (tan(screen.fov_W*(M_PI/180.0)) * camera.right[0]);
	pw[1] = screen.distance * (tan(screen.fov_W*(M_PI/180.0)) * camera.right[1]);
	pw[2] = screen.distance * (tan(screen.fov_W*(M_PI/180.0)) * camera.right[2]);
	fprintf(stdout, "(%f,%f,%f) (%f,%f,%f) \n",ph[0],ph[1],ph[2],pw[0],pw[1],pw[2]);


	for(i=0; i < screen.H; i++){
		for(j=0; j < screen.W; j++)
		{
//			CreateRay(&ray, &camera, &screen, i, j);
//			fprintf(stdout,"(%d,%d)  o (%f,%f,%f) d (%f,%f,%f)  \n ",i,j,ray.startPoint[0],ray.startPoint[1],ray.startPoint[2],ray.direction[0],ray.direction[1],ray.direction[2]);
			ds=0.0;

			dr[0] = ((screen.distance * camera.look[0]) + ph[0] + pw[0] - (i*ih));
			dr[1] = ((screen.distance * camera.look[1]) + ph[1] + pw[1] - (j*iw));
			dr[2] = ((screen.distance * camera.look[2]) + ph[2] + pw[2])	     ;
			normalizeVectoriR(dr);
//			fprintf(stdout,"(%d,%d) (%f,%f,%f) \n",i,j, dr[0], dr[1], dr[2] );	
			
			float hit = Intersection2(volume, sr, dr, &ds, &screen.map[k]);			
//			float hit = Intersection2(volume, ray.startPoint, ray.direction, &ds, &screen.map[k]);			
//			destroyRay(ray);
			//fprintf(stdout,"(%d,%d)=%f\n",i,j,ds);
			/*
			screen.map[k]	= 0;//(hit*255)/182;
			screen.map[k+1]	= hit <10.0 ? 0 : ((ds*255)/20)+100;
			screen.map[k+2]	= 0;
			*/	
//			screen.map[k]	=  ds*255;//(35<hit) 		? 255*ds : 0;
//			screen.map[k+1]	=  0;//((25<=hit)) 			? ds*128 : 0;
//			screen.map[k+2]	=  0;//((0<hit)&&(25>hit)) 		? ds*128 : 0;
			
			screen.map[k+3]	= 255;
//			if(hit)fprintf(stdout,"hit %f (%d,%d,%d)\n",hit,screen.map[k],screen.map[k+1], screen.map[k+2]);
			k+=4;
		}
			fprintf(stdout, " %d - ",(i*100)/screen.H);
	}

	toPNG(screen.map,screen.H,screen.W);

	return 0;
}
